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I. INTRODUCTION 

Modern calculations of the electronic structure of con- 
densed matter rest on density functional theory'^ (DFT), 
whose implementation is carried out within various forms 
based on the Kohn-Sham (KS) formulation of the theory 
While there is a vast body of work that attests to the utility 
of this approach, it has well documented failures in predicting 
electronic structure and related properties of materials where 
the effects of the Coulomb interaction are judged to be par- 
ticularly strong. Such materials include, but are not limited 
to, 3 d- transition metal oxides, 4/-electron rare earths and 5/- 
electron actinide systems. An important factor in this failure 
is the presence of the well-known unphysical self-interaction 
terms in the Hartree (classical) expression for the Coulomb 
energy as well as in the exchange and correlation energy func- 
tional. While there have been numerous attempts* to cor- 
rect for the presence of self-interaction, no fully satisfactory 
coherent scheme has yet emerged (for a at least a partial com- 
pendium of methods, see"^). 

Self-interaction was introduced into the theory through the 
original formulation of the so-called local density approxima- 
tion (LDA) or local spin-density approximation (LSDA)''"'', 
including their gradient corrected versions (GGAs), in which 
the Coulomb energy of the interacting particles is expressed 
as that of a classical charge distribution, n{r), with itself. 

Iff n(ri)n(r2) 
2 7 7 |ri - r2| 

An advantage of the classical form (Hartree term) is its ex- 
plicit dependence on the density that allows the functional dif- 
ferentiation with respect to the density and the determination 
of its contribution to the potential. Even though the problem- 
atic nature of the Haitree term was realized from the begin- 
ning, most of the cuiTently used energy functionals'** are still 
based on it. The presence of self-interaction reduces not only 
the reliability and convincing power of results obtained in the 



LDA, but also affects the formal standing of the theory as a 
whole. 

Among the functional forms cited in"*, some* stand out 
as strong candidates for the development of a general theory 
for a self-interaction free formulation of KS-DFT. Prominent 
among these are the self-interaction coiTection (SIC) method**, 
and the optimized effective potential (OEP) method*'^ '^. The 
former corrects for self-interaction on an orbital by orbital ba- 
sis, but cannot be shown to remove the SI error completely**. 
The latter relies on the treatment of so-called orbital depen- 
dent functionals and attempts a functional differentiation of 
the Coulomb energy by means of a chain rule based on the 
w-representability of the density. 

It is well known that self-interaction does not arise when 
the Coulomb energy is calculated in terms of the pair density, 
a quantity that in the KS formulation of DFT is determined 
through the Slater determinant of the KS orbitals. However, 
the dependence of the pair density on the density is only im- 
plicit. In one implementation of the OEP method, this im- 
plicit behavior is attacked through the chain rule of func- 
tional derivatives which, in turn, involves the calculation of 
the (inverse) susceptibility of the KS system^. Alternative 
procedures, such as the parametrization of the potential '^ '"^ 
have been formulated, where the coefficients are determined 
to minimize the energy. Results obtained in the two methods 
are compared with those of the present procedure in a latter 
section. 

In this paper, we provide a new and alternative formula- 
tion that eliminates self-interaction effects by construction, 
preventing their appearance through the treatment of the ex- 
change term that is functionally differentiated with respect to 
the density to produce the exchange potential. Combined with 
the so-called Hartree term, this leads to a Coulomb potential 
that is determined as the functional derivative of the Coulomb 
energy, properly expressed in terms of the pair density, with 
respect to the density. We refer to the new method throughout 
the paper as self-interaction free (SIF). The formalism relies 
on a novel mathematical procedure in which a function, e.g.. 
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the KS orbitals, that depends implicitly on the density, can be 
expanded in an orthonormal and complete basis of functions 
whose elements are expressed explicitly in terms of the den- 
sity. As demonstrated in the paper, functional derivatives of 
the occupied KS orbitals can be obtained through term by term 
differentiation of the expanded forms, not requiring the calcu- 
lation of the susceptibility or its inverse. The final forms of the 
functional derivatives are expressed analytically in terms of 
the gradients of the occupied orbitals, whose numerical eval- 
uation is straightforward. 

A short introduction of the method has been presented in 
previous work-". In the following pages, we provide complete 
details of the technical, that is algebraic, and computational 
components of this method that show both the simplicity of 
the formalism as well as the power that derives from it. We 
demonstrate the characteristic features and efficacy of the new 
method by means of model one-dimensional calculations as 
well as calculations of the exchange potential and energies of 
the ground states of atomic systems from Helium to Krypton. 
We compare some of our results with the corresponding ones 
obtained in the OEP and Hartree-Fock (HF) methodologies, 
and comment on the differences between the results of the 
calculations. 

This is a rather long paper because of our intention of pro- 
viding as complete and detailed an exposition of the method 
as possible. The paper takes the following form. 

In Section II, we provide a derivation of the KS equations 
identifying explicitly the terms corresponding to the correla- 
tion energy and potential, the exchange potential determined 
in a local approximation to the KS equations, and point out the 
difficulties associated with self-interaction. Subsection II C 
states the quantum mechanically correct form of the Coulomb 
energy to be used in the KS formalism. The equidensity ba- 
sis, the formal and computational foundation of the method- 
ology introduced here is set forth in Section III. The appli- 
cation of the formalism to model and realistic atomic system 
is presented in Section IV. Differences between the method 
provided in this paper and others are discussed in Section V. 
Section VI contains our conclusions. 

In the interests of completeness, we also provide in the form 
of a set of extended appendices detailed derivations of key 
expressions resulting from our formulation that we used in 
calculating the results presented in the main part of the paper. 



II. KOHN-SHAM EQUATIONS 

In order to clarify the SIF methodology and its effective- 
ness in treating a self-interaction free Coulomb energy, we 
briefly review the basics of KS density functional theory''^. 
We consider a finite number, N, of electrons interacting via 
a Coulomb repulsion confined in a volume, il, and moving 
under the action of an external potential, v{r). For simplic- 
ity, in the following we consider only non-degenerate states 
and in much of the development we suppress the presence 
of spin. Final expressions, however, are given in full spin- 
resolved form (see appendices). 



A. Kohn-Sham Equations for Ground States 

The Hamiltonian describing an interacting system of N 
electrons in an external potential takes the usual form. 



T 



N 



u 



N 



(2) 



with the operators i), and corresponding, respectively, 
to the external field, the kinetic energy and the inter-particle 
interaction (Coulomb repulsion). The ground-state energy of 
the system is given by the expectation value. 



(3) 



where IvE*^) denotes the many-particle ground state of 
. For electrons, I'J'^) leads to a wave function, 
11, r2, . . . , rjv), that is antisymmetric with respect to in- 



terchange of the coordinates (and spins) of individual parti- 
cles, according to the requirements of Fermi statistics. We use 
the notation, \^^) — > ?^(r), and say \'^^) leads to n(r), to 
denote the property. 



n{r) = N / |'J'^^(ri,r2,...,rjv)rdr2...dr 



N, 



(4) 



where ri(r) denotes the single-particle density function nor- 
malized to the total number of particles, N . This property is 
formally equivalent to taking the expectation value with re- 
spect to ^I^^ (ri , r2, . . . , rjv) of the single-particle number op- 
erator, n(r) = ?/'^(r)7/'(r), where ijj^ and ?/; are field creation 
and destruction operators for an electron at r. We now write 
(3) in the form. 



w(r)ng(r)dr 



= Min 

n(r) 



7;(r)n(r)dr + F[n] 



= Min£;[n], (5) 

n(r) 



in terms of the constrained search functional^ 



F[n] = Min (*|T 

|*)^«(r) 



N 



(6) 



Given a density, n(r), the constrained search examines all an- 
tisymmetric A^-particle wave functions that lead to the den- 
sity and delivers the state (in the absence of degeneracy) 
that produces the minimum value of (vp^jT^ + U^\^^). 
Cioslowski""* has provided a formal procedure for generating 
all antisymmetric wave functions leading to n(r) and identify- 
ing that '^^ (ri, . . . , rjv) (in the absence of degeneracy) that 
determines F[n]. For u-representable densities, gives the 
Hohenberg and Kohn functional', F^afXri], and the minimiz- 
ing state 1 5'^) coincides with |^^)- 

For any other anti-symmetric state (wave function) |$^) 7^ 
1*^) such that 1$^) ^ n{r), we have. 



F[n] < ($'^|r + C/|$ 



(7) 
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so that the exact ground-state energy. Eg, forms a lower bound 
of the expectation values of the Hamiltonian with respect to 
antisymmetric A^-particle states, |$^) — ?^(r)- 

As in the initial formulation of DFT by Kohn and Sham"", 
we postulate the existence of a fictitious non-interacting A^- 
particle system described by the Hamiltonian, 



H 



N 



(8) 



under the action of an external potential, Vg, whose ground- 
state density is identical to the density of the interacting sys- 
tem described by . In analogy with (6), we define the 
constrained search functional. 



r,[n]= Min ($^|T^|$^). 

|$")^n(r) 



(9) 



In the absence of degeneracy, the minimizing |$^) is a single 
Slater determinant^^ (denoted by the subscript s) of order N. 
Because generally, \^^) ^ 1'^^), we have. 



Es = J v{r)n{r)dr + Fs[n] > Eg, 



(10) 



where 



F,N = ($f|r^ + [/^|$0> ^W- (11) 



Defining the quantity, Ec[n] = F[n] — Fs[n], and adding and 
subtracting Fs[n] to E[n] in Eq. (5), we obtain. 



E[n] 



w(r)n(r)dr + Fs[n] + Ef.[n]. 



(12) 



The quantity Ec[n] is referred to as the correlation energy. 
We can view the last expression as a means of determining Eg 
through the states of the non-interacting system, given 
the functional difference, Ec [n] . 

The Slater determinant, 1$^), is obtained from the solu- 
tions of a single-particle Schrodinger equation, which then 
alos defines the potential, Vs{r), 



1. 



/z(r) = ej,(r). 



From (4) it follows that. 



N 



(13) 



(14) 



where the orbitals, fj{v), coiTespond to the N eigenvalues of 
(13) that lie the lowest in energy. With respect to the same 
states, we also define the pair density for the non-interacting 
system. 



1 ^ 

2El/*(^i)/^(r2)-/,(ri)/.(r2)r (16) 

i,3 



y"|$^(ri,r2,...,rw)f dr3...dr^ (15) 



Now, the functional Eg [n] takes the form. 



Mr) 



1-2 



(17) 



The form of Fs[n] is expressed in terms of the fully ex- 
changed two-particle density and is by construction free of 
self-interaction effects. From Eq. (13) we obtain the expecta- 
tion value of the kinetic energy operator. 



1. 



Mr) 



N 



N 



J dr/;(r)/,(r)- J drv,{r)n{v) 



dr Vs{r) n{r), 



(18) 



where a star denotes the complex conjugate of a quantity. 
Using these expressions, we can write 

E,ln] = E[n] - Es[n] 
= (T[n]-r,[n]) 

+ I dr, I dr, <r,,v,-[n\)-nAv,,r,-[n]) 

where T denotes the exact expectation value of the kinetic en- 
ergy of the interacting system and 7i(ri, r,) is the correspond- 
ing exact two-particle density. 

Using the stationarity property of the KS energy for the 



ground state with respect to the density 



1^ <5£;g[Ti] 



5n{r) 







and 



ST, [n] 
(5n(r) 



-Vs{v) + c^'* (where in the following the con- 
stant c is supressed), we obtain the requirement at the ground 
state. 



Vsiv) = v{v) + 



dvi / dr 



Sn, (ri,r2) 
5n{r) 

\ri - r,! 



A*c(r) 



(20) 

where the terms in the square brackets are the sum of the 
Coulomb potential and the correlation potential. 



6Er_ 



5n{v) 



(21) 



completing the expression for Vg to be used in Eq. (13). Equa- 
tion (20) gives a unique KS potential, leaving no freedom in 
its determination: It is the unique (generally within a constant) 
potential that leads to the density. 

The quantity Ec [n] appearing in the previous discussion is 
exact, but generally unknown. For the case of systems with 
infinite numbers of electrons, a correlation energy functional 
that is consistent with the second Hohenberg-Kohn theorem 
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in yielding energies that are no lower than their experimental 
counterparts can be constructed based on the properties of the 
homogeneous electron gas (jellium), as developed in the work 
of Ceperley and Alder-''. For finite systems, such a consistent 
construction of the correlation energy is still lacking. In the 
following, we set this term equal to zero and concentrate on 
what is known as the exchange only form of the KS functional 
(the nomenclature to become clear below), that takes the form, 



E\n] 



v{r)n{r)dr + Fs[n] > Eg 



(22) 



Our aim is to develop a methodology that, with /ic(r) set 
equal to zero, determines the potential (KS potential), Ws(r), 
within the exchange-only mode of implementation of the the- 
ory, through the derivative of the Coulomb energy with respect 
to the density including the exchange term. 



B. Hartree Approximation of the Coulomb Energy 

In the local-density approximation (LDA)^, the Coulomb 
energy of an interacting A^-particle system, e.g., an atom, is 
approximated by the classical expression. 



1 



n(ri)n(r2) 
|ri - rzl 



dridr2, 



(23) 



where n(r) denotes the ground-state density of the interact- 
ing system, a form that lends itself immediately to functional 
differentiation with respect to n{r), (the density appears ex- 
plicitly inside the integral). 

Computational ease, however, comes at a high price: the 
product of densities allows the simultaneous occupation of the 
same orbital by a single electron, in violation of the Pauli ex- 
clusion principle. Alternatively, an electron at position ri is 
allowed to interact with itself at r2, leading to a clearly un- 
physical self-interaction ejfect. 



C. Coulomb Exact-Exchange Kohn-Sham Functional 

Within the KS formulation, the quantum-mechanically cor- 
rect expression of the Coulomb repulsion energy of an non- 
interacting A^-particle system takes the form. 



If ns(ri,r2 

C^QM = / rdridr2, 

2 J ri-r2 



(24) 



where ris (ri, r2) is the pair density obtained from the KS or- 
bitals. In general, it is convenient to define the exchange term, 

-^s(i"i,r2) = ns(ri,r2) - n(ri)n(r2), (25) 

and a corresponding exchange energy. 



1 /• Js(ri,r2) ^ , 

o 1 rdridr2. 

2 J ri-r2 



(26) 



In this case, with fj (r) denoting an orbital in the KS deter- 
minant, the pair density takes the form (a denotes spin). 



"-slri,r2 
1 
2 



jEl/»(^i)/j(^2)-/,(ri)/.(r2)|' (27) 



fi.a, (l"l)/j.<T, i^2)Aa, (1-2) 



and the exchange term has the form, 

Js(ri,r2) = 



(28) 



E 



fL. (^l)fla, iT^2).fj.a, (ri)/.,^, (r2) (5^, 



(29) 



Finally, the exchange contribution to the single-particle po- 
tential arising from the Coulomb energy, called the exchange 
potential Vx, can be written as. 



Wx(r) 



6n{r) ' 



(30) 



in which the functional dependence on the density is explicitly 
indicated. 

Written in terms of the expression in (27), the energy func- 
tional to be minimized in the exchange-only form of the KS 
formulation of ground-state DFT takes the form, 

E[n] = J v{r)n{r)dr + T,[n] + Uci[n] + Js[n]. (31) 

This form exhibits clearly the difficulty that has been en- 
countered in attempts to implement the KS formulation of 
DFT. While the functional differentiation of the Hartree term 
with respect to the density is straightforward (explicit depen- 
dence on density under the integral), the orbitals that enter the 
definition of the exchange term are only implicitly dependent 
on the density defying immediate differentiation with respect 
to n{r) by analytic means. Even a brief survey of the many 
forms proposed for bypassing this seemingly impossible task 
would be exceedingly lengthy and take us too far afield of our 
intended purpose. At the same time, a few comments are in 
order. 

The OEP method*'^ is the best known and most widely 
used procedure for the solution of the KS equations in the 
presence of exact exchange. The functional differentiation of 
the Coulomb energy in reconstructing the potential within the 
OEP method is computationally intensive, usually involving 
the calculation of the inverse susceptibility and ultimately re- 
quiring the solution of an integral equation. Alternatively, the 
potential can be parametrized'^ '** where parameters are de- 
termined in order to minimize the energy. The present SIF 
method avoids the need to solve the OEP equations, by using 
a strictly analytic treatment in obtaining the functional deriva- 
tive of the Coulomb energy with respect to the density. This 
is accomplished through the use of the equidensity basis for- 
malism, introduced by Macke^^ and Harriman^*". In addition 
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to orthonormality, Zumbach and Maschke^^ proved that the 
basis is complete. The expansion of the orbitals in this basis 
brings out the expHcit dependence on the density allowing the 
straightforward differentiation using the product rule. The ex- 
istence of the basis has been known for some time, but as far 
as we are aware the use proposed here has not been attempted. 



III. EQUIDENSITY BASIS 

The calculation of the Coulomb potential in Eq. (20) re- 
quires the functional derivative of the non-interacting pair 
density ns(ri,r2) with respect to the (spin) density n°'(r). 
This can be performed analytically by expressing the former 
in terms of an expansion that exhibits explicitly the single- 
particle density. Namely, for each orbital, //(r), we write. 



integral sign in Eq. (29). The procedure is as follows: We re- 
place the orbitals in the exchange term with their expansions 
in terms of the equidensity basis, and use the property of func- 
tional derivatives. 



(37) 



(32) 



to perform the functional differentiation of the exchange term. 
It is clear that any density can be used to construct the 
equidensity basis, and the corresponding expansion of the 
KS orbitals. However, only the density corresponding to the 
KS orbitals is useful to construct the equidensity basis. This 
choice allows a direct functional differentiation of the basis 
with respect to the density at each iteration step. 



A. Coulomb Potential Calculation 



'^k(r, K]) 



n'^(r) 



TV" 



exp{ik-R'"(r, K])}, (33) 



where 4)^{y, [n'^]) are the elements of an orthonormal and 
complete basis^^"^^, the equidensity basis, where k — 
{kx, ky,kz) denote a set of three signed integers, the a^''^ are 
expansion coefficients, and the vector R'^(r, [n°']) is defined 
by the expressions. 



i??(x,2/,z,K]) = 



27r 



iV-(y,z,[n-]) 



2n 



oo 



(34) 



with 



N''{y.,z,[n"]) ^ I dx rf{x' ,y,z,[n'']) 

J —OO 

/OO 
dy'7V-(y',z,K]) 
-OO 
rOO 

N^ln"] = / dz' N^iz^in"]) (35) 



where < R2, R3 < 2tt. The transformation, r — > R, 
maps three-dimensional coordinate space onto the volume of a 
cube of side 27r with the points at infinity mapped onto the sur- 
face of the volume. Note that the R are explicit functionals of 
the density. This particular choice of i?i , i?2, R3 is not unique, 
with other choices, e.g. permuting coordinates or other coor- 
dinate systems, being possible^** as well. In the following we 
use the definitions given in equations (34) and (35). 
The coefficients, a^'^, are given as the overlap integrals, 

«r = / /;(r)0r(r,W) dr. (36) 

The functional derivative of the Coulomb energy can now be 
carried out through the differentiation of the orbitals under the 



In general, the Coulomb potential is given by the functional 
derivative of the Coulomb energy with respect to the density. 



Sn{r') 
1 



dri dr2 C/(ri,r2) ns{ri,r2] 



dridr2 t/(ri,r2) 



(5ns(ri,r2) 



(38) 



For three-dimensional systems the interaction U is given by 
the expression. 



f/(ri,r2) = 



1 



|ri - r2| 



(39) 



(where clearly U has no functional dependence on the den- 
sity). 

We split the pair density into two parts, the well known 
Hartree and the exchange contributions (with z in a com- 
pound index that includes spin), 



?^6(rl,r2) 



N 



E 



/;(ri)/;(r2)/.(ri)/,(r2) 



-/;(ri)/;(r2)/,(ri)/,(r2)<5..,., 
7i(ri) n{v2) 



N 

E 



/*(ri)/;(r2)/,(ri)/.(r2)<5.„.^. 



n(ri)n(r2)-l- Js(ri,r2) 



(40) 

(41) 
(42) 



where the exchange term has contributions from likewise 
spins only. In this form, the sum runs over all i and j. Since 
the potential from the Hartree contribution is trivial, in the 
following we concentrate on the functional derivative of the 
exchange term, Js(ri, r2). 
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Applying the product rule of functional differentiation yields 
the expression 

^ SJs{ri,r2) 
^ (5n'^(r') 

N 

= -E 



+ 



^/;(r2)/,(ri)/,(r2) 



+/nri)/;(r2)/,(ri)||^ 



(43) 



In general, the integrals obtained combining (38), (42) and 
(43) have the form (U is symmetric in ri and r2). 



J J dridr2C/(ri,r2)/:(r2)/,*(ri)/c(r2) 



Sfdjri) 
Sn{r') ' 



(44) 



where a, b,c,d G We define the integral : 

P'in) := I dr2C/(ri,r2)/;(r2)/,(r2), (45) 

that is known to be smooth. The detailed evaluation of 
these integrals for three-dimensional systems is given in Ap- 
pendix A. 

Since U is symmetric in the spatial coordinates and = 
[P'^y equation (43) can be further simplified. We use the fact 
that ri and r2 as well as the summation indices i and j can be 
interchanged, to obtain the expression. 



<(r')= o 



dri dri C/(ri,r2 



(5Js(ri,r2) 



N 

E 



/dr/-(r)/,(r)0« 



-'(7i,a-i .(7 



drPHv) f*(r) -^M^ 



(46) 



<(r') - - 



25R f: / dr (r) /; (r ) , (47) 



where 6c,i.aj,a equals 1 when <Ji = aj = a and vanishes 
otherwise. The exchange potential is always real. It becomes 
clear that the crucial quantity of interest is the derivative of an 
orbital with respect to the density, 

B. Explicit Derivatives of Orbitals 

Within SIF, the derivative of an orbital with respect to the 
density is straightforward. Using the expansion (32) and the 



definition of the equidensity basis (33) we obtain the expres- 
sion, 

E^l^-Ej^orcc-.i"-]) (48) 

:.R S^n-^{r) 



E 



1 



Sn''{r') 

6- 1 



(Sn'^(r') 



I til O Iv • ; 



6a^ 



(5n<^(r') 



E 



k L 



(5(r-r') i 



1 



+i a 



(r, K^])k 



+ 0kHr, |r^"' 



^n'^(^') 



5{r - r') /f* 



2n'^i(r) 



E 



2iV'^- 



(5n°'(r') 



(49) 



For the sake of formal completeness, we include the func- 
tional derivative of the normalization integrals N°' with re- 
spect to the density, that leads to the second term after the last 
equals sign. 

When the spin of the orbital and the spin of the density with 
respect to which it is differentiated coincide, the last expres- 
sion reduces to the form, 

<5/f(r) <5(r-r')^.^^^_ /f(r) 



SR{r, [n" 

5a, 



2N° 



k 



(50) 



where the second term (from the derivative of the normaliza- 
tion) leads to a constant shift in the potential and can be ne- 
glected. In the following we also neglect the last term con- 
taining the functional derivative of the expansion coefficients 
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with respect to the density. The reason for making this seem- 
ing approximation will be discussed in other work. 

The third term of Eq. (50) involves an infinite sum over k, 
raising questions about convergence. This problem is circum- 
vented through the realization that the same infinite sum oc- 
curs in the gradients of the orbitals. As shown in Appendix C, 
comparing functional and spatial derivatives allows the re- 
placement of the sum over k with expressions involving par- 
tial derivatives of the orbitals. For example, as shown there, 
in the one-dimensional case we obtain the result. 



2tt n'^{x) 
N" 1 



1 



2n<^{x) 



27r y^n^Jx) \ y^n'^ix) 



(51) 



where primes on functions denote spatial derivatives. The 
three-dimensional case follows analogously, leading to the ex- 
pression. 




2-iTN''{y,z) 

N" 
y 2ttN''(z) 



V,/ 



/(r) 


dn'iT) 


2Ti"(r) 


dx 


/(r) 




2n^{r) 


9y 


f(r) 




271" (r) 


dz 



QxPl2 



-'23 



(52) 



with Pab being the partial derivatives of R with respect to the 
coordinates. 



dRg 
db 



a, be {x,y,z}, 



(53) 



given for the three-dimensional case in Appendix B 1. 
In short, we can write 



<5/f(r) _<5(r-r') 



(5R(r, [n 



/r(r) 



Sn(r') 



(54) 



where the Q include summations over k to infinite order. 

Still to be considered is the functional derivative of R with 
respect to the density. The derivation is shown in detail in 
Appendix B 2. For the one-dimensional case we obtain the 
expression. 



5Ri{x,[n'^]) _ 27r 



Q{x-x') 



Sw^ix') N° 
while in three-dimensions we find 



(55) 



dR2{y,z,[n-]) 
dn'^ir") 



dn<^{r") 



2i:Q{x-x")- Rl{x,y,z) 

5{z - z") 
N'^{z) ^ 

27rQ{y -y")-RUy,z) 



1 



27re(z 



.(56) 



We have now derived an analytic and closed-form expression 
for the functional derivative of an orbital with respect to the 
density. Through the use of the product rule, these expressions 
can be used to obtain the contribution, Wx(r), to the Coulomb 
potential. Computational details and the final expressions for 
the three-dimensional case are given in Appendix D. 

Before closing this section it should be noted that our re- 
sults for the exchange energies and potentials neglecting of the 
last term in Eq. (50) match those of Harbola and Sahni"^, as 
discussed later in Section V C, which are obtained from elec- 
trostatics, through a method that is exact for spherical charges 
and independent of an expansion in an orthonormal basis. 



C. Formal Summary 

We summarize the discussion of the previous section. The 
use of spatial derivatives allows us to eliminate the explicit 
evaluation of the equidensity basis and the expansion coeffi- 
cients, leading to a closed-form expression for the functional 
derivative of an orbital with respect to the density. In the one- 
dimensional case we obtain the expression. 



Sn'^{x") 2n'^{x) 



d(x — X ) 



+ 



n'^{x) 



2N' 
r{x) n'^'ix) 



Q{x - x").(57) 



2n'^{x) 

The analogous results in three-dimensions take the form. 



2N° 



(58) 



[27re(x- a;") - Rl{x,y,z)\ 



5{y~y")5{z-z") 
y'^' r2^e(y-2/")-i?^(y,z)] 



N'^iz) 



5{z-z") 
[2^Q{z-z")-Rl{z)\, 



(59) 



with the quantities R defined in (34), the N in (35) and the 
Q in (52). Computational details are given in Appendix D. 
The calculation of the potential becomes now straightforward 
using equations (47) and (45). 
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D. Iterative Procedure 

The solution of the KS equations rehes on an iterative pro- 
cedure, with an updated KS potential and density determined 
at each iteration step. The sequence of steps in the treatment 
of the KS equations is summarized below: 

1. At each step, i, of the iteration, determine the orbitals, 
/j*''(r), and the density, rS^\v). 

2. Determine the derivatives of the orbitals with the re- 
spect to the density, take the functional derivative of the 
Coulomb energy with respect to the density in terms of 
spatial gradients, and obtain the Coulomb energy con- 
tribution to the KS potential. 

3. Solve the KS equation for the new potential, go back 
to the first step and iterate until convergence is reached 
within some preset tolerance. 

In short, the only difference with conventional procedures is 
the treatment of the full Coulomb potential expressed in terms 
of the pair density, rather than just the Hartree term or modi- 
fications to it. 



IV. EXAMPLES 

In this section we present the results of applications of our 
method to two systems for which the exact answers are known 
analytically, to a one dimensional model system, and finally 
on the series of atomic systems from Helium to Kiypton. For 
the latter case, we compare our results to those from HF and 
OEP methods and discuss the differences between the results 
of the three approaches. 



A. Analytic Examples 

The formalism allows the expression of the functional 
derivatives of the exchange energy with respect to the den- 
sity by analytic means and leads to closed-form expressions 
in terms of the spatial gradients of the orbital functions. These 
expressions can be compared to exact results in cases where 
analytic expressions are known. 



1. Two-Electron Systems 

As a first example we discuss the ground state of a two- 
electron systems, such as the Helium atom or the Hydrogen 
dimer (-^2), with two electrons of opposite spin in the same 
spatial orbital under the same external potential. 

Assuming real, and node free orbitals we obtain 



i=l 

fir) = y/w^ir) 



(60) 
(61) 



Sfjr) ^ 1 
Sf 1 



S{r - r') 



(62) 
(63) 



from which the exchange potential can be determined as 
follows: 

<i-) = \J /dridr2f/(ri,r2)'^:|^ (64) 



fin 



dri dr2 C/(ri,r2; 
/(r2)/(r2) 



+/(r2)'^/(ri)/(r0 



(65) 



dri dr2 [/(ri,r2; 
^<5(r-ri)/(r2)/(r2) 
-t-i5(r-r2)/(ri)/(ri) 
= - J dri C/(ri,r)n'^(ri) 

<(r) = J dri C/(ri,r)n(ri) 

= _lj/Hartree/ \ 
2 ^ 



(66) 



(67) 



The exchange potential is exactly half of the Hartree potential 
but with the opposite sign, so that the self-interaction error is 
half of the Hartree term. 

Using the method proposed in this paper, we write 

/(r) = J2 ak0k(r) with ao = 1, (68) 

k=0 

where all coefficients other than k = vanish. From Eq. (50) 
we obtain the expression, 

<5/(r) _6{r-r')^^^^_ 5{r-r') ^ ^^^^ 



Sn'^ir') 2n'^(r) ' ' ' 2/(r) ' 
Using Eq. (47) for the exchange potential yields the result. 



(70) 



<(0 = -2/dr/"(r)/(r)^ 

= -/"(r') = - y dr U{r, r')/(r) /(r) (71) 

= - J dr C/(r,r')n'"(r) (72) 

= -^/dr C/(r,r')n(r) = ^Hartree(^,)(73) 

For this simple example, the SIF method reproduces the cor- 
rect analytic expression for the exchange potential. 
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Hartree Potential 



Our second example deals with the Hartree term. The func- 
tional derivative of the Hartree energy with respect to the 
density is well known. From 



E" = 



1 



dr dr' 



, n(r)n(r') 



we obtain 



5E 



H 



5n{r") 
It can also be shown, that 

5E" 5E" 



V^«(r") 



5E 



dr- 



H 



n{v) 



(74) 



(75) 



5n{r") (5n1-(r") (5n^(r") 
We now express the density in terms of orbitals. 



l/^(r"). (76) 



E 

i=l 



/;(r)/.(r) 



(77) 



and express the orbitals in terms of the equidensity basis, and 
take the functional derivative of E^ written in terms of the ex- 
panded forms. The results leads to . The detailed deriva- 
tion of this result is shown in Appendix E. 



B. Numerical Examples 

In this section we apply our method to one-dimensional sys- 
tems in terms of the particles-in-a-box problem and to realistic 
atomic systems. 



1. One-dimensional Square Well 

The work reported in this section is designed to test a sim- 
ple case of non-interacting, spinless Fermions confined in an 
one-dimensional well of length, L = \xi — xq\ = 1, with infi- 
nite potential walls. We choose iV = 6. This example is used 
for illustrating the method, deriving the potential correspond- 
ing to the energy for a given form of the inter-particle inter- 
action. We restrict ourselves to a non-self consistent solution 
for a vanishing (or constant) potential. We choose an inter- 
particle interaction that decays exponentially with respect to 
inter-particle distance. 



(78) 



that allows us the freedom of manipulating the range of the 
interaction and assess its effects on the exchange potential. 



Eq. (14) and the pair density, Eq. (42). The normalized wave 
functions of this system are given by the expressions. 



fn{x) 



(^) 0<.<L, (79) 



with quantum numbers rt = l,2,3,.... The energies read as 



En — 



2^ virJ 



(80) 



with the ground-state density given by n{x) = 

In analogy to Eq. (34), we define the quantities 
Ri/2/3{x, y, z), that in one dimension reduce-' to a function 



q{x) = ^ / n{x')dx'. 



(81) 



For the six-electron case, q{x) is shown in Fig. 1. We see that 




FIG. 1. The quantity q(x) (solid blue line) for the one-dimensional 
square well example with 6 electrons. The red dashed line shows 
q{x) = X which corresponds to a constant density. In this case the 
basis functions become plane waves. 

the function q{x) is almost linear An exact linear behavior, 
q{x) = X, (red dashed line) would correspond to a constant 
density with the equidensity orbitals reduced to plane-waves. 
The functional derivative of f'^!'^}-, takes the form, 

on{x ) ' 



6n{x") 



2ti 
TV 



Q{x-x"). 



(82) 



We also have. 



The lowest in energy orbitals in this case are known an- 
alytically and lead to an analytic expression for the density in 
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In spite of the dependence on k, the sum on the right-hand side 
converges sufficiently rapidly to be numerically stable. In this 
case, the equidensity basis functions, (t)k{x), take the form. 



21 functions 



0fc(a;) 



"-(^) ikq{x) 

N 



n{x) 



^ikq(x,[7i]) 



(83) 



There exist two choices for k in constructing a complete and 
orthonormal set of basis functions: Either signed whole in- 
teger values, fc = 0, ±1, ±2, ±3, . . ., or half-integer values, 
k = ±i,±|,±|,.... The rate of convergence of the expan- 
sion of the orbitals is found to depend on spatial symmetry, 
with those even under reflection (symmetric) about the center 
of the box being described more efficiently by whole integer 
values of fc, while half-integer values lead to faster conver- 
gence of the orbitals that are odd under reflection (antisym- 
metric). A plot of the coefficients vs. values of k is shown in 
Fig. 2. We find that with the choice of the faster converging 




FIG. 2. Log-log plot of expansion coefficients of spatially symmet- 
ric (even under reflection) and antisymmeytric (odd under reflection) 
orbitals in terms of an equidensity basis, as indicated in the legend. 




FIG. 3. Ground-state density of six non-interacting electrons in an 
one-dimensional box with infinite walls as discussed in the text. The 
blue dashed line marks the analytic expression, plotted below the 
other curve. At this scale they are almost indistinguishable. 




expansion the error in the norm becomes smaller than 10^^ 
when more than 100 basis functions are taken into account. 
Generally, we find that fewer than 1,000 functions are suffi- 
cient for convergence. 

The ground-state density of the system, (the six electrons 
occupying the orbitals labeled n — 1, 2, . . . , 6), is shown in 
Fig 3, while the pair density is shown in Fig. 4, with the num- 
ber of basis functions indicated in the panel. Within the res- 
olution of the figure, the two results are essentially indistin- 
guishable. The same rate of convergence characterizes the 
pair density, shown in Fig. 4. We note the vanishing of the pair 
density along the line xi — X2, as expected from Eq. (40). 



FIG. 4. Non-interacting Ground-state pair density of six non- 
interacting electrons in an one-dimensional box with infinite walls 
as discussed in the text. 



Another quantity to look at for convergence is the infinite 
sum over fc appearing in Eq. (50). The result with the sums 
carried out to infinite order is given by the analytic expres- 
sion in Eq. (51). The sums over fc for the first three orbitals 
are shown in Fig. 5 for different numbers of expansion co- 
efficients and basis functions taken into account, with those 
used chosen symmetrically around fc = 0. For this example. 
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results shown here correspond to 1,000 terms taken into ac- 
count. The analytic expressions lead to the same result. An 
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FIG. 5. Infinite sum, see Eq. (51), evaluated for the first three orbitals 
of the six electron example. Convergence is shown by increasing 
the numbers of coefficients (chosen symmetrically around 0) in 
comparison to the analytic result. 



we find a rather good agreement using a few tens of coeffi- 
cients in a broad range in space, between about a; = 0.1 and 
X = 0.9 in Fig. 5. At the boundaries, the difference between 
the analytic expression and the truncated sums increases, the 
oscillations being a manifestation of the Gibbs phenomenon. 
Even though the results show convergence with the number 
of terms taken into account, we use the closed-form expres- 
sions of Eq. (51). This bypasses the explicit construction of 
the equidensity basis and the expansion coefficients. 

Fig. 6 shows the potential obtained for the four values of 
A = 1, 10, 50, 100, as defined in (78), using both the Haitree 
expression (blue dashed curves) for the interaction energy, 
Eq. (23), as well as the quantum mechanically correct expres- 
sion, Eq. (24) (red solid line). In all cases, convergence is 
well established with about 200 basis functions, although the 





FIG. 6. Single-particle potential contributed by the interaction en- 
ergy in the system of six non-interacting electrons confined in a box 
with infinite walls as discussed in the text. The blue dashed line is 
calculated from Eq. (23) (Hartree), whereas the red solid line is based 
on Eq. (24) (pair density). 

immediate effect of the quantum expression are the lower val- 
ues of the potential obtained from Eq. (24), caused by the pres- 
ence of the exchange term that reduces the effective region of 
interaction of two electrons. Of interest is also the behavior of 
the two different expressions for the potential in the limit of 
short-range interaction, or large values of A. The Hartree term 
yields a potential that approaches the form of the density, es- 
sentially reaching that form for A — 100 (upper curve in panel 
on lower right-hand corner). In that limit, the interaction po- 
tential resembles a delta-function in two-particle space and the 
energy given by the Hartree term is simply proportional to the 
density at the point xi = X2- 

The behavior is considerably different when the pair density 
is used to calculate the interaction energy. Now, in the limit of 
short-range interaction interaction, and spinless fermions, two 
paiticles cease interacting altogether because the pair density 
vanishes in the region in which the interaction is noticeably 
non-zero and the contribution to the single-particle potential 
drops essentially to zero (lower curve in panel in lower right- 
hand comer in the figure). 



2. Atomic Calculations 

In the following we apply the SIF method to realistic 
atomic systems. Details of the calculations, in particular of 
the evaluation of the Coulomb and exchange potential are 
given in Appendix D. Because the form of the correlation en- 
ergy is cuiTently unknown, the results presented are all fully 
self-consistent and converged within the exchange-only mode 
(Ec = 0). This allows us to compare SIF results with those 
obtained within the exchange-only OEP and the HF methods. 

Fig. 7 show the results for Helium, Lithium, Beryllium, 
Boron, Neon and Argon obtained from a self-consistent 
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exchange-only calculation. The left column shows the ex- 
change potential on a logarithmic scale that allows a detailed 
study of the behavior close to the nucleus. To emphasize the 
asymptotic behavior, in the middle panel, the same exchange 
potential is shown on a linear scale, but multiplied by the dis- 
tance to the nucleus, r. Finally, the density is plotted in the 
right column. For all examples, the corresponding full OEP 
results are shown for comparison. Lithium and Boron are 
spin-polarized, resulting in different potentials for the major- 
ity and minority spin channel. 

For Helium, the exchange potential within SIF coincides 
with the one obtained by the OEP method, and therefore the 
densities are also identical. This identity is not preserved for 
heavier atoms. As seen in Fig. 7, the SIF potential has always 
a non-negative slope, whereas the OEP potentials exhibits re- 
gions where the slope becomes negative. In the OEP litera- 
ture^ this is referred to as the intershell structure. In Section V, 
we provide a formal explanation of this behavior. 

Even though the potentials are different, the SIF and OEP 
densities are indistinguishable at the scale of the plot. All ex- 
change potentials behave regularly near the origin, and have 
an — 1/r asymptotic behavior, as expected. For illustrative 
purposes, we plot the exchange potentials for the atoms He- 
lium to Kiypton in Fig. 8. 

Total energies for the atom series obtained in SIF, OEP and 
HE are compared with experiment in Table I and are plotted 
as difference to experiment in Fig. 9. Experimental data for 
the total energy can be obtained by summing all ionization po- 
tentials of the corresponding atom, a process whose difficulty 
increases prohibitively with atomic number, Z. Published re- 
sults for the total energy are often characterized by a number 
of corrections to account for zero-point motion and relativistic 
effects. The resulting rather unclear situation for the total en- 
ergies is discussed in Appendix F. In any case, the experimen- 
tal values of the ground-state energies lie the lowest compared 
with those of HE, OEP and SIF in corresponding increasing 
order. We discuss this order in more detail in the following 
section. 



V. DIFFERENCES BETWEEN HE, OEP AND THE SIF 
METHODS 

A. Comparison with Hartree-Fock 

It has become customary in the literature to compare the 
results of an exchange-only implementation of DFT to those 
of the HE method, and we follow this practice. 

The HE method relies on the direct calculation of a determi- 
nantal wave function to minimize the energy of an interacting 
system. The single-particle orbitals defining the determinant 
are the basic variables of the method, and their unrestricted 
variation provides a path to the lowest energy obtained in this 
procedure. The method, however, does not provide a path to 
the determination of the ground state of an interacting system. 

By contrast, the basic variable in DFT is the density, and 
the KS orbitals are restricted to those that lead to the density 
in terms of a non-interacting collection of electrons in an ef- 



fective potential, Vs{r). In order to minimize the energy the 
functional derivative with respect to the density is taken, and 
evaluated at that particular density formed by the KS orbitals. 
In an exchange-only implementation of the theory, the restric- 
tion of the orbitals to reproduce the density is expected to lead 
to higher values of the total energy compared with those of 
the HE method. This, indeed, is the case in all calculations 
reported here, see Table 1. The restrictions can also be dis- 
cussed in terms of the potential. Within the KS scheme, the 
potential, see Eq. (20), is local, as in HE the potential becomes 
orbital-dependent and non-local. 



B. Comparison with OEP 

Since the most widely used method for obtaining a poten- 
tial from an energy functional written explicitly in terms of 
the orbitals is the OEP, we compare our results with the ones 
obtained in that method. Unfortunately, published results, e.g. 
for atomic systems, obtained by the OEP method are not con- 
clusive. There are at least two different approaches in imple- 
mentations of the OEP. The first one is based on the original 
work by Talman and Shadwick''', where the susceptibility is 
approximated by first order perturbation theory and then in- 
verted. 

A different approach was suggested by Wu and Yang'^ '*^ 
where the potential is expanded in a basis and the expansion 
coefficients are varied in attempting to reach convergence. In 
principle, both approaches should determine the OEP poten- 
tial and energies uniquely. The uniqueness of the OEP results 
has been discussed in more detail in the literature''^''''*. Numer- 
ically, however, the two approaches seem to differ as indicated 
in Fig. 10 showing the exchange potential for Ar determined 
by four different methods. Heaton-Burgess, Bulat, and Yang''^ 
examined in great detail the numerical convergence of their 
method and concluded that the oscillatory behavior in the ex- 
change potential (around r = 10^^ and 1 in the curve marked 
OEP), referred to as the inter-shell structure, is unphysical. 
The same interpretation of the inter-shell structure as unphys- 
ical has been proposed by Harbola and Sahni (see below). 

In their study, Heaton-Burgess et al.''^, found that the ampli- 
tude of oscillations is related to the finite basis set, or in other 
words, the oscillations are the result of an unconverged treat- 
ment. For Argon, their converged exchange potential is sim- 
ilar to SIF, and does not show any of the so-called inter-shell 
structure, characterizing the results based on the conventional 
treatment'^, where the functional derivative of an orbital with 
respect to the potential is treated within first-order perturba- 
tion theory. 

The relative order of the energies obtained in HE, OEP and 
the SIF methodologies can be understood in terms of the vari- 
ational constraints characterizing each of the methods. As 
mentioned above, the unrestricted HE method (in the absence 
of a correlation-energy functional) yields the lowest energies, 
while the variation of the wave function (KS orbitals) re- 
stricted to reproduce the density yields the highest. The OEP 
method employs a local potential (hence gives energies higher 
than HE ) but allows variations in terms of the orbitals treated 
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z 


Symbol 


Element 


Exp(") 


EXp('')(Ref.'") 


x-SlF 


x-OEP(Ref."') 


HF(Ref.") 


2 


He 


Helium 


-2.90339 


-2.904 


-2.862 


-2.862 


-2.862 


3 


Li 


Lithium 


-7.47798 


-7.47806 


-7.432 


-7.433 


-7.433 


4 


Be 


Beryllium 


-14.6685 


-14.66736 


-14.571 


-14.572 


-14.573 


5 


B 


Boron 


-24.6582 


-24.65391 


-24.527 


-24.528 


-24.529 


6 


C 


Carbon 


-37.8557±0.0002 


-37.8450 


-37.687 


-37.689 


-37.690 


7 


N 


Nitrogen 


-54.6119 


-54.5892 


-54.401 


-54.403 


-54.405 


8 


O 


Oxygen 


-75.1100±0.0001 


-75.0673 


-74.809 


-74.812 


-74.814 


9 


F 


Fluorine 


-99.8062±0.0001 


-99.7339 


-99.406 


-99.409 


-99.411 


10 


Ne 


Neon 


-129.0507±0.0059 


-128.9376 


-128.542 


-128.545 


-128.547 


11 


Na 


Sodium 


-162.4283±0.0032 


-162.2546 


-161.852 


-161.857 


-161.859 


12 


Mg 


Magnesium 


-200.3 100±0.0033 


-200.053 


-199.606 


-199.612 


-199.615 


13 


Al 


Aluminum 


-242.7 121±0.0037 


-242.346 


-241.868 


-241.873 


-241.877 


14 


Si 


Silicon 


-289.8683±0.0037 


-289.359 


-288.845 


-288.851 


-288.854 


15 


P 


Phosphorus 


-34L9464±0.0087 


-341.259 


-340.709 


-340.715 


-340.719 


16 


S 


Sulfur 


-399.0351zhO.0050 


-398.110 


-397.495 


-397.502 


-397.506 


17 


CI 


Chlorine 


-461.3813±0.0051 


-460.148 


-459.470 


-459.478 


-459.483 


18 


Ar 


Argon 


-529.1122±0.0093 


-527.540 


-526.804 


-526.812 


-526.817 


19 


K 


Potassium 


-601.9677±0.0515 




-599.150 


-599.159 


-599.165 


20 


Ca 


Calcium 


-680.1022±0.0679 




-676.743 


-676.752 


-676.758 


21 


Sc 


Scandium 






-759.718 


-759.728 


-759.736 


22 


Ti 


Titanium 






-848.375 


-848.397 


-848.407 


23 


V 


Vanadium 






-942.852 


-942.876 


-942.886 


24 


Cr 


Chromium 






-1043.334 


-1043.350 


-1043.360 


25 


Mn 


Manganese 






-1149.848 


-1149.860 


-1149.870 


26 


Fe 


Iron 






-1262.425 


-1262.440 


-1262.450 


27 


Co 


Cobalt 






-1381.376 


-1381.410 


-1381.420 


28 


Ni 


Nickel 






-1506.828 


-1506.860 


-1506.870 


29 


Cu 


Copper 






-1638.938 


-1638.950 


-1638.960 


30 


Zn 


Zinc 






-1777.820 


-1777.830 


-1777.850 


31 


Ga 


Gallium 






-1923.235 


-1923.250 


-1923.260 


32 


Ge 


Germanium 






-2075.335 


-2075.350 


-2075.360 


33 


As 


Arsenic 






-2234.215 


-2234.230 


-2234.240 


34 


Se 


Selenium 






-2399.844 


-2399.860 


-2399.870 


35 


Br 


Bromine 






-2572.416 


-2572.430 


-2572.440 


36 


Kr 


Krypton 






-2752.029 


-2752.040 


-2752.050 



TABLE 1. Total energies in Hartree for the atom series from Helium to Krypton. Two columns with experimental data are given, EXP^"' 
andEXP**"', for details see Appendix F. Other columns show the total energies obtained by various methods (SIF, OEP and HF ) within the 
exchange-only mode. The data are plotted in Fig. 9. In general, the SIF energies for the atom series lie slightly above the OEP ones, and both 
of them are higher than the HF results. The differences are discussed in the text. 



as independent variables leading to energies lower than SIF C. Alternative Method 

where the orbitals are required to reproduce the density. This 
feature is further exemplified through the use of unconverged 

(and finite) basis set that leads to OEP energies similar to those For special cases, e.g., closed shell atomic systems and gen- 
ofFLF "*^. erally spherical charge densities, Harbola and Sahni^** have 

shown that the exchange potential can be obtained through 
the use of classical electrostatics. Notably, their method does 
not require the calculation of the susceptibility or its inverse 
and is independent of both SIF and OEP. For the case of Ar, 
the exchange potential found using their procedure precisely 
matches that which we find using the SIF approach (see Fig- 
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FIG. 7. Exchange potentials for Helium, Lithium, Beryllium, Boron, Neon and Argon. The left columns depicts the radial exchange potential 
on a logarithmic scale, the middle column the exchange potential multiplied by r, the distance from the nucleus, so as to exhibit its asymptotic 
behavior, and the third column the density. Results of this paper (referred as x-SIF) are compared to OEP calculations. The curves represent 
self-consistent results obtained in the absence of a correlation energy (exchange only). 
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FIG. 9. Total energy differences, normalized by the number of elec- 
trons, between SIF, x-OEP, and HF to the experimental values, see 
Table I. a) refers to EXP'°\ and b) to the column EXP''''. For com- 
parison, LSDA energies'" are shown in panel a). 



FIG. 8. SIF exchange potentials multiplied by r, the distance from 
the nucleus, for the atomic series Helium {Z — 2) to Krypton (Z = 
36). Both figures show the same data, just at a different angle. The 
colored surface is for illustration, and exchange potentials are only 
meaningful for integer Z values, marked by the blaclc lines. 



ure 10). It is also noteworthy that Harbola and Sahni provide 
a physical interpretation of the negative slope of Vs obtained 
in some OEP calculations. It signifies the unphysical property 
that a negative charge, an electron, is repelled by a positive 
exchange hole. 

In later work, Sahni^^^ ""* states that the OEP results for the 
exchange potential are not the functional derivative of the ex- 
change energy, but include an additional contribution to the 
energy, to which he refers as correlation-kinetic effects. 

In still further work'''', Harbola and Sahni determined the 
ground-state energies of atomic systems from He to Rn. For 
closed shell systems, where densities are spherical, their re- 
sults, that are exact in this case, coincide with our results 
shown in Table I and Fig. 10. 



VI. CONCLUSIONS 

The developments in the paper sprung from a particularly 
simple observation: The functional derivative of the Hartree 
term with respect to the density is facilitated by the explicit 
appearance of the density in the integrand. The same ease 
would materialize were the orbitals in the exchange term given 
explicitly in terms of the density in a differentiable functional 
form. The expansion of each orbital in the equidensity basis 
is particularly convenient for satisfying this requirement. 

Based on this expansion, we have proposed a novel treat- 
ment for the Coulomb energy to be used in implementations 
of the KS formulation of DFT within a local approximation. 
In this treatment, the Coulomb energy is expressed in terms of 
the pair density constructed from the single-particle orbitals 
arising from the solution of the KS equations, thus avoiding 
by construction self-interaction effects. The pair density, in 
turn, is written in closed-form explicitly in terms of the den- 
sity. This is accomplished through an expansion of the single- 
paiticle orbitals in terms of the equidensity basis whose ele- 
ments are written explicitly in terms of the density. We have 
demonstrated how infinite sums appearing in the formalism 
can be replaced by quantities containing spatial derivatives 
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Appendix A: Calculation of the integrals 

In spherical coodinates, the integrals defined in Eq. (45) 
can be evaluated using the well known formula based on an 
expansion of U into spherical harmonics, 

' ' L={l,m) > 



FIG. 10. Comparison between the exchange potential of Argon atom 
derived by the methods used in this paper (SIF, red, solid), the ap- 
proach of Harbola and Sahni"' (dashed, blue), and the results ob- 
tained by two different implementations of the OEP method, as indi- 
cated: OEP designates the results obtained within first order pertur- 



and 



bation theory (dashed, green), and Heaton-Burgess et al.' 
orange) those of directly varying the potential. 



(dashed. 



and can thus evaluate the sums (Eqn. (51) and (52)) to infinite 
order. The procedure is consistent with the bounds of DFT 
guaranteeing that the energies calculated within the method 
will be no lower than those of the true ground-state energy. In 
contrast to other methodologies directed at the determination 
of the exchange potential, such as the OEP, we find that the ap- 
proach taken here is consistent with the fundamental premise 
of DFT: the only independent variable is the density, and it is 
the only variable with respect to which functional differentia- 
tion in principle can be performed. 

We have shown that the SIF method provides a computa- 
tional simple formalism to obtain the exchange potential as an 
upper bound of the energy, avoiding the calculation of the sus- 
ceptibility and its inverse or the solution of an integral equa- 
tion, where here only KS orbitals and their spatial derivatives 
are needed. Results for atomic systems are very close to those 
obtained by the OEP method and match the ones by Harbola 
and Sahni. In the form discussed in this paper, the formal- 
ism can immediately be applied to all non-periodic systems, 
e.g. molecules. An extension to periodic systems (solids) is 
currently the subject of ongoing work. 



(A2) 



In particular for a spherically symmetric function g we get, 
g(r') ^, _ 471 



rdr' 



g (r' ) r'Mr' + 47r / g{r')r'dr' 



(A3) 

If the orbitals are of the form / = RL{r) Yl{6, cf)) the angu- 
lar part can be integrated analytically using the Gaunt coeffi- 
cients. 



L'L" 



dfYL'{f)Yl„{f)YLir) L^{l,m). 



Then the integrals have the form, 

P^{v,)^ £■ i?^^(ri)>L(0i>i) 



(A4) 



(A5) 



with Imax = 2max(/i, Ij) and 
(ri) = i?*(ri)i?,(ri) (c^^l)*. (A6) 



where, depending on the symmetry, some coefficients R^l 
might vanish. Since the functional derivatives of R contain 
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6 distributions it will be difficult to expand them into a spher- 
ical harmonics basis. We find it more useful to evaluate the 
integrals in Cartesian coordinates and perform the remain- 
ing integrals to obtain the potential. 



Appendix B: Deriatives of the transformed coordinates R 

1. Spatial Deriatives 

For completeness, we provide all spatial derivatives of the 
functionals R(ri,(r, [n]) defined in Eqs. (34). 



/322 := 



27r 



dR2 
dy 

dz ~ N{z) 



'R2 



N{y,z) 
N{z) 



(B5) 



27r / dy' / Ax'—n{x',y',z) 



00 



00 



dy' / dx' —n{x' ,y' ,z) 



/3ii 
P12 



dRi 

dx 

dRi 



= 2n 



N{y,z 
1 



/332 := 



(Bl) 



dy N{y,z) 

-Ri 
dRi 1 



2n 



dx'---n{x',y, z) 
oy 



P33 ■■= - 
We see that the Jacobian is 



dR3 
dx 

dR3 
dy 

dz 






27r 



d 

dx'— n(a;',y, z) 
dy 



(B2) 



det /3 = det 



dz N{y,z) 
-Ri 



27r 



d 

dx' — n{x',y,z) 
dz 



dRi dRi dRi 
dx dy dz 
Q dR2 dR2 



dy dz 







dRs 
dz 



N 



N 



(B6) 

(B7) 
(B8) 
(B9) 



n(r), (BIO) 



dx'—n{x\y,z) 
dz 



^ j^' n{r) where d is the dimension of 



dx 



(B3) 



(B4) 



and in general det [3 
the system. 

2. Functional Derivatives of the R 



Straightforward functional differentiations lead to the ex- 
pressions; 



6Ri{x,y,z, [n'^]) 
6n{r") 



= 2tt 



2tt 



^"(r") /_!lda;'n(a;',y,z,[n'^]) 
j:^dx'5{x'-x")Siy-ynSiz-z") 
j"^^dx'n(x',y,z, [n'^]) 

/OO 
dx'6{x' - x'')S{y - y")d{z ~ z") 
-00 

J^^dx' n{x',y,z, [n'^]) 
J^^dx' n{x',y,z,[n'^]) 
d{y-y")S{z-z") 



7V(y,z, K 



2tt Q{x — x") — Ri{x, y, z, n[a]) 



SR2{y,z, Kl) 
Sn{r") 



= 2tt 



2tt 

-2tt 



S /!^dy'/_°°^dx'n(x',y',z,K]) 
Snir") dy' J^^ dx' n{x' , y', z, [n'^]) 

/-co ^y' Too " ^nsjy' - y")5{z - z") 

/^dy'/^da;'n(a;',y',z, [n'']) 



dy' / dx'5{x' - x")5{y' - y")5{z - z") 

) J —00 

/-ooV/_°ldx'n(x',y',z, K]) 



!°°oo d^;' ni^x', y', z, [n'^]) 
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and 



(5n(r") 



5{z - z") 



2TTQ{y-y'')-R2{y,z,n[<j]) 



27r 



2tt 



6 11^ dz' dy' da:' nix', y',z', [n^]) 



-oo 

5n{r") dz' JZ^dy' dx' n(x', y', z' , 

/-oo d^' Coo dy' - y")s(z' - z") 



d^' dy' dx' n{x', y', z', [n-]) 



-2tt 



dz' 



dy' / dx'(5(a;' - x")S{y' - y")S{z' - z") 



/-oo d-' /-°1 dy' Coo d-' "(-', y', -', Kl) 



1 



/-^dz'/_!^d?;'/^da;'n(a;',y',z', [n'^]) 
2t: Q{z ~ z") - Ri{z,n[a]) 



Appendix C: Spatial Derivatives and their Connection to 
Infinite Sums 



The use of the equidensity basis requires the evaluation of 
the rather difficult infinite sum in Eq. (50) whose value de- 
pends on the number of terms taken into account. The follow- 
ing considerations allow the sum to be carried out to infinite 
order. 

First, we discuss the one-dimensional case. We are inter- 
ested in evaluating the sum. 



k k 

(CI) 

where the orbitals f{x) are expressed in the form (where for 
simplicity of notation we omit the spin superscript on the or- 
bitals). 



ak 



''^"i^) ikq{x,\n'']) 



From the definition. 



q{x, K]) - 



2lT 



we have. 



and 



_ 27r 



n''{x')dx' 



e{x-x'), 



Mx^ = 

dx ^ ' 



(C2) 



(C3) 



(C4) 



(C5) 



The spatial derivative is given by 

„ df{x) 



/' 



da; 



E 

k 

E 



dx 



^''■'^(x) ikq(x, In"]) 



(C6) 



--o-k 



N" 2^/vMx) dx 



+ ik ak 



n'^jx) ikqix^ln"]) dq{x) 



NO 



dx 



2n-{x)^'"'y N- 



+ ^n''{x)J2'^kik 



N 



from which it follows that 



n'^{x) 



,ikq{x,ln'']} 



(C7) 
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fix) 
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2TTn'^{x) 
N'^ 1 
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y/n^ix) 



f(x)n° 

fix) \ 
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fix) 



27r .Jn'^{x) \^/v°\x) 



(C8) 
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(CIO) 



The spatial derivatives with respect to x (marked by primes) 
of the orbitals and the density can be calculated numerically. 

It now follows that it is no longer necessary to con- 
struct explicitly the equidensity basis, or the expansion coef- 
ficients, flk. The three-dimensional generalization of the one- 
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dimensional results is straightforward. The components of the 
gradient of an orbital, /, decomposed into the equidensity ba- 
sis are given by 



and 



1 dn'^ir) 



fir) 



J k-R(rJri'']) 



(Cll) 



2n'^(r) da 

+«2^0k(r, [n ]) [ ^ k j , (C12) 



with a = {x, y, z}. 

From Eqs. (C12), (Bl), (B4) and (B7) we obtain 
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and 



Q: = E«(r)(ifci 



N^iy,z) 
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2n'^{r) dx 

In the same way, using Eqs. (C12), (B2), (B5) and (B8) 
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(C15) 



we obtain the expression. 
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Finally using Eqs. (C12), (B3), (B6) and (B9), we find, 
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^^•^ 271- (r) 9z y,/^23 
These results can be summarized in vector form. 
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Appendix D: Numerical Evaluation of the Exchange Potential 
Components 

We split the exchange potential into different components 
with regard to Eqn. (54) and (47) 

<(r') = <[ll(r') + <[Q^l(r') + v:^^y\r') + <[Q^I(r'). 

(Dl) 

In the following we show how to evaluate the different in- 
tegrals to obtain the exchange potential, in cartesian coordi- 
nates, which makes it easy to deal with integrals over Theta 
functions. 

The contribution from the first term of Eq. (54) takes the 
form, 

^[1] 



vf^{r') = -2^Y.^-^-.-. 



y'drir^"(ri)/;(ri)/,(ri) 



'^(ri - r') 
2n-(ri) 



The contribution that includes f Qx) in the last term of 



Eq. (54), reads as follows 



d n 



-2?ft£S.,,.,..a [ dr/^J"(r)/;(r) 

2ti Q{x - x') - y, z, [n"]) 

/,(r) dn-{v) 



5{y-v')5{z~z') 
N"{y,z,[n]) 



27rn'^(r) 



2n<^{r) dx 



= -2di dxdydz 



2TTe{x - x') - Ri{x,y,z, K]) 
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^ Sa, ,a {x, y, z) f* {x , y, z) 



6{y-y')5{z^z') 



Qlix,y,z) 

N^{y,z) 



(D2) 



The integrals over y and z are trivial. Left over are two simple 
integral over x. 
We can define: 



Ax(x,y,z) ^6a^.a^,a PHx,y,z) f*{x,y,z] 

ij 

Al{x,y,z) ^6a,.aj,a PHx,y,z) f*{x,y,z] 



A.l{x,y,z) := ^S^^.a^a {x,y, z) f* {x,y, z) 
that leads to 

<[Q'^i(r') = <[Q'^i(-',y',^') 

-23* / dx2TTA'^{x,y',z') 

J x' 
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N^{y,z) 
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(D3) 



The contribution from the term that includes (^|^ Qy) in 
the last term of equation (54), reads as follows 
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The contribution from the term that includes Qz) in 
the last term of equation (54), reads as follows 
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this is a constant, spin dependent 



(D6) 



Appendix E: Functional Derivative of the Hartree term 



For the Hartree term it is easy to show that 

SUh 5Uh 5Uh 



Sn^{r) 5n-^{r) dn{r) 



Vnir). (El) 



Nevertheless, it still needs to be proven that the proposed pro- 
cedure calculating the functional derivative with respect to the 
density leads to the exact same result, namely the Hartree po- 
tential. 

For the spin density we have 



i 



(E2) 



(E3) 



where a E {x, y, z} and the sum runs over all orbitals of spin 
a. In the spin polarized case there is no functional dependency 
of an orbital with spin cr and the spin density with the opposite 
spin. 



(E4) 



Now we can show that the following expression vanishes (here 
shown for cr 
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showing that there is no contribution to the potential from the 
term appearing in the functional derivative of an orbital with 
respect to the density ^^^py- • "k * ^ (^k [n(r)] 
_ SRln{r)] Q npj^g remaining four terms used to calculate 
the potential contribution have the same structure, 

J J dridr2t/(ri,r2) 

ij 

x/;(r2) /.(ri) f,{r,) 5{r, - r')<5..,. (E8) 

i j 

= \\j Av2U{y',V2)n{v2) ^\v"(v'). (E9) 

Due to the product rule during the functional differentiation 
this term appears four times so that the contribution to the 
potential is exactly the Hartree potential. One should note 
that taking the functional derivative with respect to the spin- 
up or spin-down density leads to the same result. All the above 
relationships apply in the same way in the non-spin polarized 
case. 

We have now proven that the proposed procedure to calcu- 
late the functional derivative of the Hartree term leads exactly 
to the Hartree potential. 



Appendix F: Experimental Data for the Atom Series 



In practice, it is not easy and in most cases impossible to 
measure quantities corresponding to the total energy of cal- 
culations. In principle, one has to ionize all electrons and 
measure the energy. The total energy is then the sum of all 
ionization energies. For the first few atoms of the periodic 
table this seems to be possible, but with larger Z it becomes 
harder and harder to ionize all of the electrons. A complete set 
of wavelength corresponding to all of the ionization potentials 
from Hydrogen to Calcium {Z ^ 20) has been published'"' 
quite some time ago. One can now use the CODATA'*'"'*^, 
containing the latest values of physical constants and their er- 
rors, to convert the measured wave length to energies. The 
second column of Table 1 labeled EXP'"' contains these data, 
with the estimated error bars. 

Based on the experimental data'"' some effort has been 
made to try to calculate non-relativistic correlation energies 
and relativistic corrections to the ionization potentials'"* "^^ 
The data from reference^" seems to be widely used in the 
literature. But since they were determined by calculations 
rather than measurements, or only partially, we have doubts 
that these energies can be quoted as experimental data. Be- 
cause it is not clear which data set should be used, we decided 
to quote both in this paper. 
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